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Abstract 

Recent experimental and theoretical studies suggest that crystallization and glass-like solidifica- 
tion are useful analogies for understanding cell ordering in confluent biological tissues. It remains 
unexplored how cellular ordering contributes to pattern formation during morphogenesis. With 
a computational model we show that a system of elongated, cohering biological cells can get dy- 
namically arrested in a network pattern. Our model provides a new explanation for the formation 
of cellular networks in culture systems that exclude intercellular interaction via chemotaxis or 
mechanical traction. 
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By aligning locally with one another, cells of elongated shape form ordered, crystalline 
configurations in cell cultures of, e.g. fibroblasts [TJ |2], mesenchymal stem cells [2J, and 
endothelial cells [3]. Initially the cells form small clusters of aligned cells; the clusters then 
grow and the range over which cells align increases with time (21 H] . To study the emergence 
of such crystalline cellular ordering, it is useful to make an analogy with liquid crystals 
[2j. For example, a "cellular temperature" can be defined to describe the cell-type specific 
persistence (low cellular temperature) or randomness (high cellular temperature) of cell 
motility, where cells of high cellular temperature (e.g., fibroblasts) are less likely to form 
crystalline configurations than cells of low temperature (e.g., mesenchymal stem cells) [2J. It 
was similarly proposed that collective cell motion in crowded cell sheets can be understood as 
system approaching a glass transition [51 E] . Although these studies provide useful insights 
into the ordering of cells in confluent cell layers, it remains unexplored how crystallization 
and glass-like dynamics contribute to the formation of more complex shapes and patterns 
during biological morphogenesis. 

Cells' organizing into network-like structures, as it occurs for example during blood vessel 
development, is a suitable system to study how cellular ordering participates in pattern 
formation. In cell cultures after stimulation by growth factors (VEGFs, FGFs), endothelial 
cells elongate and form vascular-like network structures [THH] . The mechanisms that drive 
the aggregation of endothelial cells and their subsequent organization into network is a 
subject of debate. Most models assume an attractive force between cells, either due to 
chemotaxis |T0H18J or due to mechanical traction via the extracellular matrix [TU1- 24j . In 
vitro experiments show that astroglia-related rat C6 cells and muscle-related C212 cells can 
form network-like structures on a rigid culture substrate [25], which excludes formation of 
mechanical or chemical attraction between cells. Therefore a second class of explanations 
proposed that cells form networks by adhering better to locally elongated configurations of 
cells [25] or elongated cells [2H]. Here we show that, in absence of mechanical or chemical 
fields such mechanisms are unnecessary: elongated cells organize into network structures if 
they move and rotate randomly, and adhere to adjacent cells. As the cells align locally with 
one another, a network pattern appears.. Additional cell-cell attraction mechanisms, e.g., 
chemotaxis or mechanotaxis, act to stabilize the pattern and fix its wave length. 

To model the collective movement of elongated cells, we use the cellular Potts method 
(CPM) [23 [2H], a lattice-based, Monte-Carlo model that has been used to model devel- 



2 



Round 
Chemotaxis 



Long 
Chemotaxis 



Long 



A 




v 



FIG. 1: Effect of chemotaxis and cell shape on pattern formation. A round, chemotacting, 

and adhesive cells (10,000 MCS), B elongated, chemotacting and adhesive cells (10,000 
MCS), and C elongated, non-chemotacting and adhesive cells (250,000 MCS). In all panels 
700 cells are seeded on the center 500x500 pixels of an 800x800 lattice. 

opmental mechanisms including somitogenesis [29], [30], convergent extension [31] and fruit 
fly retinal patterning [32]. The CPM represents cells as connected patches of lattice sites 
with identical spin a G N + ; lattice sites with spin a = represent the extracellular matrix 
(ECM). To simulate stochastic cell motility, the CPM iteratively displaces cell-cell and cell- 
ECM boundaries by attempting to copy the spin of a randomly selected site into a randomly 
selected adjacent lattice site x, monitoring the resulting change AH of a Hamiltonian, 



A copy attempt will always be accepted if AH < 0, if AH > a copy attempt is accepted 
with the Boltzmann probability P(AH) = exp(AH / ^(a)) , with /i(cr) a "cellular temper- 
ature" to simulate cell-autonomous random motility. For simplicity, we here assume that 
all cells have identical temperature. The time unit is a Monte Carlo step (MCS), which 
corresponds with as many copy attempts as there are lattice sites. 

The first term of Eq. [T] defines an adhesion energy, with the Kronecker delta returning a 
value of 1 for site pairs at cell-cell and cell-ECM interfaces, or zero otherwise. In the model 
two contact energies are defined: J ce ii,ceii for a > at both lattice sites, and J ce ii,ECM for 
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a = at one lattice site. The second and third term are shape constraints that penalize 
deviations from a target shape, with A and L a target area and length, and a(a) and 1(a) 
the current area and length of the cell; Aa and Al are shape parameters. We efficiently 
estimate 1(a) by keeping track of a cellular inertia tensor as previously described |14j . 

In a subset of simulations, we further assume that cells secrete a diffusing chemoattractant 
c, which we describe with a partial differential equation: 

dc(x, t) 



DV 2 c(x, t) + s(l - 5(a(x), 0)) - e 5(a(x), 0)), (2) 



with diffusion constant D, secretion rate s and decay rate e. After each MCS, a forward 
Euler method solves Eq. [2] for 15 steps with At = 2 s with zero boundary conditions. To 
model the cells' chemotaxis up concentration gradients of the chemoattractant, during each 
copy attempt from x to x' we increase AH with a Aif chemotaxis = A c (c(x) — c(x')), with 
A c a chemotactic strength [33J. One lattice unit (l.u.) corresponds with 2/im. We use the 
following parameter settings, unless specified otherwise: // = 1; J C eii,ceii = -5; Jceii,ECM = -35; 
Aa = 1; A L = .1; A c = 10; A = 100 l.u. 2 ; L = 60 l.u.; D = 10~ 13 mV 1 ; e = 1.8 ■ 10" 4 s~ x ; 
s = 1.8 • 10 -4 s -1 . Unless stated otherwise, a simulation is initialized with 175 cells randomly 
distributed on a 220x220 area at the center of a 400x400 lattice. 

As Fig. [l] shows, and in agreement with previous reports [H], if we allow for chemotaxis, 
rounded cells accumulate into rounded clusters (Fig. [l]A.) and elongated cells aggregate into 
networks (Fig.[ljEJ). Interestingly, however, chemotaxis is not required for network formation: 
cell-cell adhesion between elongated cells suffices for forming networks (Fig. [l)3). Movies 
corresponding with Fig. [IjEJ and C [38] suggest that the gradual alignment of cells with 
their neighbors is key to network formation and network evolution. To characterize this 
cell alignment, we define 6(x,r) as the angle between the direction of the long axis v(a(x)) 
of the cell at x, and a local director n(x), a weighted local average of cell orientations 
defined at radius r around x: n(x,r) = (v(a(y))) ^ €Z2 .^_^ <r y- Figure [2]A. and B depict 
the value of 6(x, 3) for simulations without (A) and with chemotaxis (B), with dark gray 
values indicating values of 6(x, 3) — > 7r/2. Network branches are separated by large values of 
6(x, 3), indicating that within branches cells are aligned, whereas branch points are "lattice 
defects" in which cells with different orientations meet. 

Supplemental Movies S3 and S4 [39] show how the cells align gradually over time in 
the absence and presence of chemotaxis. To characterize the temporal development of cell 
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FIG. 2: Crystalline cell ordering during network formation. A-B 6(x,r) with r = 3 for a 
simulation with (A) and without chemotaxis (B) after 25,000 MCS. C Temporal evolution 
of orientational order parameter 5(r) for r = 20 (black curves), r = 40 (gray curves) and 
r — > oo (light gray) without chemotaxis (solid) and with chemotaxis (dashed). Order 
parameter is averaged over 10 simulation repeats (gray shadows represent standard 

deviation) . 



alignment in more detail, we use an orientational order parameter 5(r) = ^ 3cos gMgVlrri 
[34] with c(cr) the center of mass of cell a. 5 ranges from for randomly distributed cells 
to 1 for ordered cells. 

Figure |2]C shows the evolution of the global orientational order parameter lim^oo 5(r) 
and of the local orientational order parameters 5(20) and 5(40). Both with chemotaxis 
(dashed lines) and without (solid lines), 5(20) grows more quickly and reaches higher order- 
ing than 5(40). The reason for this is that in cells of length 50 — 60 l.u., 5(20) (covering cells 
up to a radius r = 20 from the cell's center of mass) only detects lateral alignment of cells, 
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whereas a radius 5(40) also detects linear line-up of cells. Thus cell-cell adhesion of long 
cells quickly aligns cells with the left and right neighbors, while it aligns them more slowly 
with those in front and behind. This results in networks with short branches of aligned cells. 
Interestingly, chemotaxis aligns cells more rapidly, both along the short and long sides of 
cells, resulting in networks with much longer branches than with adhesion alone. 
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FIG. 3: Relation between cluster size and cell displacement. Clusters are calculated for 
each morphology between 500 and 25,0000 MCS (100 simulation repeats), with an interval 
of 500 MCS; see text for details. The error bars represent the standard error of the linear 

fits used to estimate diffusion coefficients. 

Next we analyze the mechanisms that drive the orientational ordering in the cell net- 
works. Visual inspection of the simulation movies suggests that single cells move and rotate 
much more rapidly than locally aligned clusters of cells. A network of locally aligned cells 
forms rapidly from initially dispersed cells. Merging of branches seems to be a much slower 
process, and potentially prevents a further evolution to global nematic order. To quantify 
these observations we measured the translational and rotational diffusion coefficients of cells 
for clusters of a given size. The translational diffusion coefficient, D t , derives from the mean 
square displacement (MSD) of a set of cells: MSD(i) = ((c(a, t) — c(a,0)) 2 ) a = AD t t. Simi- 
larly, the rotational diffusion coefficient, D T , derives from the mean square rotation (MSR) 
of a set of cells: MSR(t) = ((a(a,t) — a(cr,0)) 2 ) a = 2D r t, with a(a, t) — «(cr, 0) the angular 
displacement of a cell between time and t. We loosely define a cluster as a connected set 
of at least two cells with relative orientations < 5°, i.e., in Fig. |2]A and B dark gray values 
separate the clusters. To detect clusters computationally, we first identify the connected 
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sets for which 6{x, 3) < 5°, which are surrounded by lattice sites of a = or sites with 
9(x, 3) > 5°. We then eliminate connected sets of fewer than fifty lattice sites. The CPM 
cells sharing at at least 50% of their lattice sites with one of the remaining sets form a clus- 
ter. During the evolution of a simulation, cells may move between clusters, and clusters can 
merge. Therefore, to calculate the MSR and MSD of cells as a function of cluster size, for 
100 simulations of 250,000 MCS we calculated trajectories of each individual cell and kept 
track of the size of the cluster it belonged to. We then split up the trajectories to obtain 
cell motility parameters for bins of cluster size of five cells wide, and used those to calculate 
D t and D r . 

The translational diffusion, D t , increases slightly with cluster size (Fig. [3]A.). This may 
reflect that the probably of hopping between small clusters will be larger than the probability 
of hopping between larger clusters, resulting in an overrepresentation of slow cells in the 
small clusters. Interestingly, the rotational diffusion D r drops with the cluster size (Fig. 
|3b), indicating that cells in large clusters rotate slower. These results suggest that the 
rotation of cells in big clusters is limited, which reduces the probability that two clusters 
rotate and merge into a single larger cluster. Therefore, if the size of clusters increases, 
their rotation speeds drop as does the probability of cluster fusion. Thus, although further 
alignment of clusters would reduce the pattern energy H (Eq. [I]), the pattern evolution 
essentially freezes. 

To corroborate our hypothesis that network patterns are transient patterns that increas- 
ingly slowly evolve towards nematic order, we looked for model parameters that could speed 
up pattern evolution. Fig. [I]A. shows the effect of surface tension (7 ce ii,ECM) on the ability 
of cells to form networks after 100,000 MCS, as expressed by the configuration's compact- 
ness C = 4 efiUl , where ^4huii is the area of the convex hull of the largest connected group 

-^hull 

of cells, and A ce \\ s is the summed area of the cells inside the hull. A value of C — > 1 
indicates a spheroid of cells, where for networks C would tend to zero. For values of 
7ceii,ECM = </ C eii,ECM — Jccl ^ cc " > o, the equilibrium pattern should minimize its surface area 
with the ECM. Indeed at increased surface tensions the cells settle down in spheroids or 
networks with only few meshes, although they initially still form network-like patterns (see 
Movie S5 [3D]). To confirm that also for 7 ce ii,ECM = 0.1 (i.e., the values used in Figs. 1- 
3) spheroids are stable configurations, we initialized our model with a spheroid (Fig. [4^3). 
Although initially some cells sprout (Fig. WC) from the spheroid due to their elongation, 
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they then align gradually and the cell cluster remains spherical. No network formation was 
detected in simulations of 100,000 MCS (Fig. [4^3) , suggesting that spheroids represent the 
global minimum of the Hamiltonian. Interestingly, in presence of chemotaxis networks form 
for a wide range of surface tensions (inset Fig. [4] and |14j). 

Our analysis suggests that in the cellular Potts model elongated, adhesive cells can form 
networks in a parameter regime where a spheroid pattern is the minimal energy state. The 
cells initially align with nearby cells, thus forming the branches of the network. In order 
for the pattern to evolve further towards the minimal-energy spheroid pattern, the locally 
aligned clusters of cells must join adjacent branches, for which they must move and rotate. 
Our analysis of the rotational and translational diffusion of cells in Fig. [3] shows that this 
becomes more difficult for cells belonging to larger clusters. Thus the networks evolve 
ever more slowly to the minimal energy state, and gets dynamically arrested in a network- 
like configuration, a phenomenon reminiscent of the glass transition, as e.g. observed in 
attractive colloid systems |33j, colloid rod suspensions [36] . and collective cell migration of 
biological cells in vitro [5J E] • Fig. [4] suggests that the cellular Potts simulations undergo 
a glass transition as the surface tension drops: for high surface tension the system evolves 
towards equilibrium, for lower surface tensions the system becomes jammed in a network-like 
state. Thus our model provides a new explanation for the formation of vascular networks 
in absence of chemical or mechanical intercellular attraction [25J. Interestingly, intercellular 
attraction via chemotaxis stabilizes the formation of networks in our simulations, and can 
drive sprouting from spheroids (not shown) . This suggests that networks are an equilibrium 
pattern of our system in presence of intercellular attraction. Nevertheless the present analysis 
of arrested dynamics provides new insight into the system with intercellular attraction: 
chemotaxis reinforces local ordering over a distance proportional to the diffusion length of 
the chemoattractant producing networks of a scale independent of surface tension |14j . 
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FIG. 4: A Relation between compactness and surface tension, with and without 
chemotaxis. The compactness was calculated at 100,000 MCS and averaged over 100 
simulations (error bars represent standard deviation). Simulations were initialized with 350 
cells on 260x260 path on the center of a 420x420 lattice. B-D evolution of a simulation 
initialized with 128 cell blob (On the center of a 420x420 grid. 
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